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Global Monte Carlo for fermions using ordered statistics 

A. Borigi 

'^University of Edinburgh, Department of Physics and Astronomy 
Mayfield Rd, JCMB, Edinburgh EH9 3JZ 

In this talk I discuss a new possibility for stochastic representation of the fermion determinant. The method 
can be used for global Monte Carlo fermion algorithms and is tested in the case of the Schwinger model. 
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1. INTRODUCTION 

The standard method for stochastic represen- 
tation of the fermion determinant is the pseud- 
ofermion method 1 . It allows to trade the com- 
putation of the determinant in favour of extra 
noise in the path integral of lattice QCD. The 
other alternative is to use the stochastic Taylor 
expansion of This suffers from the so called 
probability sign problem which is created by the 
noise [3]. The problem is less apparent if the 
method is applied to a fractional power of the de- 
terminant, a fact which is used by the Kentucky 
Monte Carlo algorithm 4 . 

In this talk I show that the sign problem cre- 
ated by the noise can be eliminated altogether if 
one uses a certain statistic from a sample of un- 
biased noisy estimators of the effective fermion 
action. As expected, a statistics which is an 
unbiased estimation of the effective fermion ac- 
tion cannot yield an unbiased estimation of the 
fermion determinant. The opposite logic suggests 
that a biased estimation of the former could lead 
in principle to an unbiased estimation of the latter 
if this is done " correctly" . This is indeed the case 
if one chooses a certain order statistic from the 
sample. I show below how this can be realized 
using the central limit theorem and the asymp- 
totic distribution of a central order statistic. 



2. USING THE CENTRAL LIMIT THE- 
OREM 

Let Xi , X2 , . . • , X„ be unbiased noisy esti- 
mators of the fermion effective action Scs = 



-Tr log(D^L>), where D is a generic Dirac. A 
particular way (and popular) to do it is using real 
(there is nothing special, one can use complex 
numbers as well) random vectors Zi, Z2, ■ ■ ■ , Zn 
such that: 



X,^~z}[logiD^D)]Z, 



(1) 



Computation of the right hand side can be done 
for example using Lanczos methods as described 
in detail in j5l6j . It is easy to show that Xi is an 
unbiased estimator of SeS- Another useful result 
(that so far went unappreciated and wasn't used) 
is the normality of the Xi distribution. This can 
be easily shown for Z2 noise vectors using the 
central limit theorem. Note that the whole idea 
relies on this result. 

3. BASICS OF ORDER STATISTICS 

Everyone has encountered at least once the 
notion of median as an unbiased estimator of 
the population mean. It is identified as the 
middle point after the sample of real numbers 
Xi , X2 , . . • , Xn drawn from the same continuous 
(to fix the idea) distribution f{x),x € R has been 
ordered in the ascending order: 



X(l),X^2), 



,x. 



(«) 



(2) 



This is possible if the sample size n is odd. In 
case n is even the median is defined to be any 
number in the interval [A"(„/2), A"(n/2+i)]- 

The median is a statistic as it is the arithmetic 
average. It is nothing particular to the median 
as an instance of a statistic calculated from the 
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ordered data. In general, X(^i^^,k = is 
called the fc'th ordered statistic. 

Distribution of can be found using a sim- 
ple heuristic. First, let fk{x)dx be the probability 
that X(^i:) is in an infinitesimal interval dx about 
X. This means that one of the sample variables is 
in this interval, k — 1 are less than x and n — k 
are greater than x. The number of ways of choos- 
ing these variables is the multinomial coefficient 
C(n; k — 1, l,rt — k). If F{x) is the cumulative 
distribution function of the sample variables then 
one gets: 

h{x)dx=[F{x)f-^[l~F{x)T-^f{x)dx (3) 

This result can be formally proven. As a starting 
classic is the book of H.A. David 7 . A recent 
introduction to the subject is provided in 0. 

4. ASYMPTOTIC THEORY 

The distribution of an order statistic can be of 
little use for direct analytical calculations of mo- 
ments for most of known parent distributions. Of 
course one can use numerical methods to compute 
numerical values to the desired accuracy. How- 
ever, as the sample size grows the order statis- 
tic distribution approaches some limiting distri- 
butions. The situation is analogous to the cen- 
tral limit theorem. In general the distribution 
is not necessarily asymptotically normal as one 
could have expected although this will be the case 
for the central order statistics that will be dealt 
with below. 

For < Of < 1, let « = \na\ + 1, where [na] 
represents the integer part of na. Let also F be 
absolutely continuous with / positive at F^^(a) 
and continuous at that point. Then the following 
result applies: 

Theorem (asymptotic distribution of a central 
order statistic). As n — s- oo, 



VH/[F-i(a)]^« 



F-\a) 



iV(0,l) (4) 



The proof can be found eg, in [5] . 

Example. Suppose that the distribution is sym- 
metric around the population mean /x. Let also 
assume that the variance cr^ is finite and /(/i) is 



finite and positive. For simplicity n is taken to 
be odd. Then the median A^(„+i)/2 is an unbi- 
ased estimator of /i and asymptotically normal. 
Further, Var[X(„+i)/2] « [Anf{^lf]-\ 

5. UNBIASED ESTIMATION OF THE 
FERMION DETERMINANT 

Let the Xi^i = l,...,n be a normally dis- 
tributed sample of unbiased estimators of the 
fermion effective action. I am interested in the 
expected value of e~^('' with i = [na] -\- 1 and a 
as above. Let X ~ N{^,a). It is easy to show 
that: 



(5) 



Hence, an unbiased estimation of the effective ac- 
tion gives a biased estimation of the fermion de- 
terminant. Now let Xa ■= {F^^{a) — ij)/cf. Then 
is asymptotically normal with mean /i-t-XaCr 

and variance 27rcr^a(l — a)e^° /n. For the ordered 
statistic one gets: 



< e 



(6) 



To have an unbiased estimator of the fermion 
determinant in terms of an ordered statistic one 
should have < e~"^(*' >= e^^. This can be 
achieved if: 



Xa — TTaa(l ~ a)e^'' /n 



(7) 



The crucial question here is: given n and a is 
there any solution of the above equality in terms 
of Xa^ To answer this question I express first a 
in terms of x^. Since {X — ii)/a ^ N{0, 1) then I 
have: 

a = P(X<x-„) = i + lerf(^) (8) 
Hence, the equation to be solved is: 

TTcr „, Xo 



(9) 



Using Taylor expansion around Xa = one gets 
na 



Xa ^ —[1 + Oixi)] 



(10) 
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For small the 0(a;^) can be neglected one gets 
a unique solution: 



Thus, using order statistics it is possible to have 
an unbiased estimator of the fermion determinant 
from a biased high estimator of the fermion effec- 
tive action. 

However, the estimator is only asymptotically 
unbiased meaning that the sample size has to be 
large enough. In practice, at least from the exam- 
ples below the asymptotic regime is reached for 
n = 40. In many application to lattice fermions 
cr is large anyway. Thus, requiring n >> a will 
satisfy the conditions of the asymptotic normality 
as well. 

Note also that a is unknown and has to be es- 
timated from the sample. This is a source of sys- 
tematic error that may bias the stochastic esti- 
mation of the determinant. To compute the error 
one has to substitute a with its estimator S in the 
formula for Xa- Then the average error is given 
by: 

< e-<'-+^ >=< e5^(-^) > (13) 

The error has to be computed in order to control 

the systematic errors. This can be done using 
for a large samples to guarantee small errors and 
then substitute in the above formula. 

6. APPLICATION TO GLOBAL MONTE 
CARLO ALGORITHMS 

The application of choice of the stochastic de- 
terminant is in Monte Carlo algorithms. I have 
tested here in case of lattice Schwinger model. To 
simulate the model I made new gauge field pro- 
posals using the gauge action which consists of 
10000 sweeps of local Metropolis steps followed 
by a global accept/reject step to include two de- 
generate flavors of fermions. Figure 1 shows the 
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Figure 1. Two flavors of Swchinger model on a 
16x16 lattice at /3 = 5 and bare quark Kogut- 
Susskind fermion mass m = 0.01 

square Wilson loops as a function of the linear 
size of the loop using three algorithms: 1) com- 
puting exactly the determinant ratios, 2) using 
the stochastic method with n = 30 noisy estima- 
tors and 3) the same as in 2) but with n = 40 
noisy estimators. This example shows that al- 
ready n = 40 estimators are enough to reach the 
exact result. 
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Xa = -T-, « 1 (11) 

4n 

In this case the sample size should be much 
greater then a and i can be calculated to give: 
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